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Abstract 

We develop, analyze, and test a sparse tensor product phase space Galerkin discretization framework for the 
stationary monochromatic radiative transfer problem with scattering. The mathematical model describes the 
transport of radiation on a phase space of the Cartesian product of a typically three-dimensional physical domain and 
two-dimensional angular domain. Known solution methods such as the discrete ordinates method and a spherical 
harmonics method are derived from the presented Galerkin framework. We construct sparse versions of these 
well-established methods from the framework and prove that these sparse tensor discretizations break the "curse of 
dimensionality": essentially (up to logarithmic factors in the total number of degrees of freedom) the solution 
complexity increases only as in a problem posed in the physical domain alone, while asymptotic convergence orders 
in terms of the discretization parameters remain essentially equal to those of a full tensor phase space Galerkin 
discretization. Algorithmically we compute the sparse tensor approximations by the combination technique. In 
numerical experiments on 2 + 1 and 3 + 2 dimensional phase spaces we demonstrate that the advantages of sparse 
tensorization can be leveraged in applications. 
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Introduction 

In this paper, we consider the numerical solution of the 
radiative transfer problem (RTP). This transport problem 
is stated on the phase space Q, = D x S as the Carte- 
sian product of a bounded physical domain D C M d , 
where d = 2, 3, and the unit ds -sphere as the parameter 
domain S of dimension ds = d — 1 = 1,2. The RTP (see 
e.g. Modest 2003) is then given as the task of finding the 
unknown radiative intensity u : Q — > R, a real function 
over the phase space satisfying 

s ■ V x u(x,s) + (k (x) + o(x)) 

u{x,s) = ic(x)Ib(x) + o{x) I <t>(s,s')u(x,s')ds', (la) 
Js 

u(x,s) = g(x,s), x e 3D, s ■ n(x) < 0 . (lb) 

We refer to Eq. (la) as the stationary monochro- 
matic radiative transfer equation (RTE), while Eq. (lb) 
constitutes inflow boundary conditions. A ray of light of 
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direction 5 is attenuated by absorption and scattering with 
the medium. In (la), k > 0 is the absorption coefficient, 
a > 0 the scattering coefficient, and <t> > 0 the scat- 
tering kernel or scattering phase function. The scattering 
phase function is normalized to j s 4>(s, s') ds' = 1 for 
each direction s. Sources inside the domain D are modeled 
by the blackbody intensity I\, > 0, radiation from sources 
outside of the domain or from its enclosings is prescribed 
by the boundary data g > 0. The vector n(x) denotes the 
outward unit normal vector which is defined in (almost 
every) point x on the boundary 3D of the physical domain. 

Due to the high dimension of the phase space, the 
nonlocality of the scattering operator, and the hyperbolic 
nature of the PDE, the efficient numerical simulation of 
radiative transfer is a challenging computational task even 
today. Still, radiative transfer as such or as a means of 
energy transfer among others is of interest in many appli- 
cations, e.g. in the fields of heat transfer (Modest 2003), 
neutron transport (Hebert 2010), atmospheric sciences 
(Evans 1998), medical imaging (e.g. Peng et al. 2011), or 
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other areas where transported particles interact with a 
background medium, but only negligibly with each other. 

In this paper, we extend the range of sparse tensor prod- 
uct discretization methods for the RTP investigated before 
(Grella and Schwab 2011a; Grella and Schwab 2011b; 
Grella 2013; Widmer et al. 2008) by a new phase space 
Galerkin framework. 

Apart from Monte Carlo methods for raytracing, the 
most popular deterministic approaches for the radiative 
transfer problem are the discrete ordinates method and 
the spherical harmonics approximation. We quote a brief 
overview of their properties from (Grella 2013). 

In the discrete ordinate method (DOM) or Sm- 
approximation, the angular domain is discretized by a 
number of fixed directions, which are inserted into Eq. (1) 
so that a system of spatial PDEs results. Without scatter- 
ing the equations for single directions are independent 
of each other, with scattering, however, they are coupled 
through the scattering integral. After the straightforward 
discretization of the angular domain, the spatial PDEs are 
typically solved using finite differences, finite elements, or 
finite volume methods. 

The DOM is popular as it is simple to implement, offers 
straightforward parallelization, and can capture directed 
radiation relatively well as some of the ordinates can 
usually be chosen freely. 

On the downside, the method can suffer from so-called 
"ray effects" (Lathrop 1968): due to the point evaluation 
in the angular domain, the scalar flux or incident radi- 
ation from small isotropic sources may appear star-like 
with rays emanating from the source into the chosen 
angular directions (Stone 2007, p. 2 and following). These 
effects occur especially pronounced in settings with low 
scattering and absorption, i.e. in optically thin media. 

An example for truncated series expansion is the 
method of spherical harmonics orPA/-approximation. The 
solution of Eq. (1) is replaced by a series of spherical 
harmonics up to some order N with spatially dependent 
coefficients. Due to orthogonality relations, the scattering 
part often decouples or couples only few terms depending 
on the scattering kernel. However, the system of PDEs for 
the spatial coefficient functions is always coupled by the 
transport part s ■ W x u. 

As low order series expansions in spherical harmonics 
do not permit a very localized resolution of the angular 
variable, the method performs best when the solutions are 
nearly isotropic in angle, which is the case in diffusive, 
so-called "optically thick" media. Then, rather low order 
spherical harmonics approximations might suffice for a 
good approximation. Indeed, the Pi method can be for- 
mulated as a diffusion equation for the incident radiation 
(Modest 2003, Sec. 15.4). For smooth solutions, the spher- 
ical harmonics method exhibits spectral convergence in 
angle (Grella and Schwab 2011a). 



On the other hand, beam-like solutions require a high 
spectral order to be resolved appropriately, leading to 
high computational complexity. In general, higher spec- 
tral orders also lead to a sharp increase in computational 
complexity when boundary conditions are to be satisfied 
(Modest and Yang 2008). 

When combined with a standard finite element or finite 
volume discretization in the physical domain D, the deter- 
ministic, numerical Sjv- and Pm -approximations exhibit 
the so-called "curse of dimensionality": the error (typically 
the L 2 -error of the solution) with respect to the numbers 
of degrees of freedom (DoF) and M s on the compo- 
nent domains D and S scales with the dimension d and d$ 
of the application problem as O with 

constants s and t. 

The first sparse finite element approximation method 
was proposed in (Zenger 1991) for the solution of Laplace 
equation in the unit square and cube. In this paper, Zenger 
developed the (direct) sparse grid approximation method 
which alleviates this curse of dimensionality: the compu- 
tational complexity is reduced, up to logarithmic terms, to 
that of a one-dimensional problem. 

The idea of sparse tensorization of finite element and 
finite difference methods was generalized by (Bungartz 
and Griebel 2004; Hegland 2003; Garcke 2007), and oth- 
ers, for the numerical solution of PDEs as well as for 
other applications where standard numerical methods are 
obstructed by the curse of dimensionality. 

Sparse tensor methods were first applied to radiative 
transfer by (Widmer et al. 2008). In that paper, the authors 
formulated a least squares phase space Galerkin sparse 
tensor approximation with hierarchical finite elements 
as discretization of the physical domain and wavelets in 
the angular domain. They proved that sufficient regu- 
larity of the solution provided, their method breaks the 
curse of dimensionality: the problem complexity reduces 
to log-linear in the number of degrees of freedom, while 
convergence rates deteriorate only by a logarithmic factor. 
However, the discretization of the scattering operator had 
not been addressed in that work. 

In earlier work (Grella and Schwab 2011a), we showed 
that the sparse tensor product method of (Widmer et al. 
2008) can also be combined with a spectral discretization 
involving spherical harmonics, resulting in a sparse Pm- 
method which also treats scattering. Boundary conditions 
were satisfied in a strong sense by introducing piecewise 
spectral functions in angle. 

Secondly we presented a sparse tensor version of the 
DOM as a sparse collocation method with a Galerkin 
ansatz in the physical domain and strong enforcement 
of the boundary conditions, while not yet accounting for 
scattering (Grella and Schwab 2011b). This sparse ten- 
sor Sjv-method was realized computationally with the 
sparse grid combination technique (Griebel et al. 1992) to 
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construct a sparse approximation to the radiative transfer 
solution. 

The sparse DOM was subsequently reformulated as a 
phase space Galerkin method with quadrature in angle 
(Grella 2013) in order to treat sparse Pn- and sparse 
S^-method in a more uniform manner. In this reformu- 
lation, we employed streamline upwind Petrov Galerkin 
(SUPG) stabilization and weak satisfacion of boundary 
conditions. Sparse S^-methods were derived as a direct 
sparse tensor method and implemented algorithmically 
via the combination technique. 

In the present paper, we derive a sparse P^- and sparse 
SA/-method from the same phase space Galerkin frame- 
work with transport stabilization and scattering. Bound- 
ary conditions are satisfied in a weak sense. In doing so 
we close a gap in the list of conceivable combinations 
of formulations regarding stabilization and type of angu- 
lar approximation. In contrast to the previous approach 
(Grella 2013), we stabilize the formulation in a different 
way and the analytical focus will be on the direct sparse 
approach. With transport stabilization and direct sparse 
approach we follow (Widmer et al. 2008) more closely, 
extending their work by treatment of scattering and weak 
satisfaction of the boundary conditions. 

Similar savings in computational effort are realized with 
other variational formulations, such as Petrov-Galerkin 
saddle point formulations (see e.g. Dahmen et al. (2012) 
and the references there). 

The outline of this paper is as follows. In Section 
'Phase space Galerkin method' we formulate the phase 
space Galerkin framework in operator form and outline 
how Pm and 5A/-methods can be derived from it. Then 
we develop full tensor and sparse tensor discretizations 
based on the framework and analyze and compare their 
convergence properties. 

Section 'Numerical experiments' presents several basic 
numerical experiments designed with the purpose of vali- 
dating and illustrating the theoretical convergence results. 

Finally we conclude this work in Section 'Conclusion' by 
summarizing and reviewing the results. 

Phase space Galerkin method 

We begin by introducing the radiative transfer problem in 
operator form. Using this compact notation we then state 
the variational formulation of our phase space Galerkin 
method and proceed to discretizations of the method. 

Operator formulation 

Problem (1) reads in operator form: Find the intensity 
a(*,s):Dx5->I such that 

Aw =/, u\ 3 n_ = g- (2) 

In this, 9 £2_ represents the inflow part of the bound- 
ary d£2 = 3D x S of the computational domain or 



phase space Q, = D x S. The inflow boundary is defined 
by 

3Q_ := {(x,s) eQ:xe T_(s)} (3) 

with the physical part of the inflow boundary 

r_(s) := {x € 3D: s- n{x) < 0}. (4) 

Correspondingly we define the physical part of the out- 
flow boundary as 

Y+(s):={x&dD:s-n{x)>0}. (5) 

The radiative transfer operator A = T + Q consists of 
the transport operator T, 

Tu := (s • V x + k)u, (6) 

and the scattering operator Q, 

Qu := erQiw := a (Id — E)w := a(x)u(x,s) 

f (7) 
— a(x) J $>(s,s')u(x,s r ) ds' . 

Here, Qi = Id — £ is the unity scattering operator, and 
E is the scattering integral operator, the integral of <t> and 
u. The source function f contains the sources of radiation 
in the domain, 

/ := Kl b , (8) 

andg is the incoming radiation on the boundary 3Q_, as 
in Sec. 'Introduction'. 

Properties of the scattering operator 

Aside from the positivity and normalization requirements 
already mentioned in Sec. 'Introduction', we assume an 
isotropic medium, i.e. <I> does not depend on x. As <& 
models the type of scattering, this assumption can safely 
be made for most applications (cf. Modest 2003, p. 268). 
Variations in the strength of scattering due to e.g. varying 
spatial density of the medium are modeled by the scat- 
tering coefficient er. As long as the following properties 
hold for almost every x, the complexity and convergence 
analysis later on could also be conducted without this 
assumption. 

Furthermore, if spherical scatterers are assumed, the 
scattering phase function does not vary with the 
azimuthal angle so that 3> only depends on the inner 
product of s and s\ From this it follows immediately that 
d>(s,s') = O(s-s') = <*>(*',«)• 

From here on, we shall take <i> to be forward 
dominant (cf. Kanschat 2008, Def. 1) if 3> («,«') = 
^^ 0 «£cos(£arccos(s • s')) with all a^ > 0. Then, one 
can show that E is positive semi-definite (Kanschat 2008, 
Lemmata 2 and 3), i.e. 

(v, Ev) L 2 (5) > 0 VveL 2 (S). (9) 
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Normalization and symmetry of <E> with respect to its 
arguments leads to normalization of the operator norm 
1 1 £ | \l 2 (S)^l 2 (S) = 1 (Kanschat 2008, Lemma 5). 

From these properties and a Hilbert-Schmidt theorem 
for integral operators (e. g. Knapp 2005, Thm. 2.4), one can 
derive that the spectrum of Qi lies in [0, 1] with an isolated 
eigenvalue An = 0, from which the next largest eigenvalue 
ki differs by a positive constant (Avila et al. 2011, Sec. 2.2). 

With the previous considerations, one arrives at the 
following properties of Q: 

Lemma 1. For any u e L?(Q.), the scattering operator Q 
as defined by Eq. (7) satisfies (cf. Avila et al. 2011, Eq. (11)) 



L 2 (Q.) 



IIQ«llz, 2 (fi) < IMIl^) \\u\\ L 2 (Q) , 



(u,Qu) L 2 (Q) > HQwII^) > °> 



(10) 
(11) 



in which the projector P 1 - maps u(x, •) € L 2 (S) to (ker Q)- 1 , 
the space orthogonal to the kernel ofQ, and X\ e (0, 1] is 
the smallest nonzero eigenvalue o/Qi. 

For a proof of (11) we refer to (Grella 2013). 
Variational formulation 

Our variational formulation will be based on a Galerkin 
finite element framework over the phase space Q. with 
stabilization applied to the operator RTP (2). 

A generic stabilized phase space variational formulation 

To begin with, we define the Hilbert space 

V:= [u eL 2 {Q.) : s- W x u el 2 (12)] (12) 
with the usual L 2 (Q.) inner product 

(u,v)l2(0) := j j u(x,s)v(x,s)dxds, (13) 

and the triple bar norm 

IIMII 2 := l|v||2 J(Q) + ||*.V,v||2 J(Q) + ||Q 1 v||2 J(Q)J veV. 

(14) 

For the weak enforcement of boundary conditions, we 
define the boundary form 

b(u, v) := (v, |« • n\u) L 2r d n_) = I I \s-n\uvdxds, 

(15) 

in which we have omitted the dependence of the outward 
unit normal n on the position x. This boundary form was 
introduced by Manteuffel et al. (2000, Eq. (2.16)). It is well- 
defined for functions v e L 2 (Q.) with finite inflow norm 



Combining (14) and (16) yields the new norm 

IMIi:=(|||v||| 2 + ||v||l) 1/2 , (17) 

which gives rise to the closed, linear subspace of V, 

Vi:={veV:||v||i<oo} (18) 

which, with the inner product related to || v|| v is a Hilbert 
space. For functions u, v e Vi, we define the bilinear form 



a(u, v) := (Rv, Au) L 2^ n) + 2b(u,v), 



(19) 



where R is a stabilization operator on the test function side 
yet to be specified. Together with the linear form 



l(v) := (Rv,f) L2(n) +2b(g,v), 



(20) 



the bilinear form constitutes the following variational 
problem: Find u € Vi such that 



a(u, v) = l(v) Vv € Vi. 



(21) 



Different ways of stabilization are conceivable and have 
been used in the literature, e. g. the least squares approach 
by (Manteuffel et al. 2000), or SUPG introduced by 
(Brooks and Hughes 1982). For our purposes here, we will 
choose the T-stabilized formulation (Grella and Schwab 
2011a) to avoid mesh-dependent quantities and the square 
of the scattering operator. More precisely, we set R = 
eT with a stabilization parameter £ that depends on the 
absorption coefficient k. 

Properties of the variational formulation 

At this point, we introduce the anisotropic or mixed 
Sobolev spaces H s,t (£2) = H S (D) <g> H'(S) as 



// s ' f (£2) := {v e L 2 (£2) : Df D£v € L 2 {Q,), 

0 < |a| < s, 0 < < ij 

with the corresponding mixed Sobolev norms 
given by 



0<jaj<s0<|/J|<i 



2 

L 2 (S2) ' 



(22) 



(23) 



|M|_ :=b(v,v) 1 ' 2 . 



(16) 



a 

Here, D« D"v denotes the weak derivative of v : DxS — > 
K of order \a \ w. r. t. * € D and order | ft \ w. r. t. s € S, with 
the multi-indices a e and /) 6 Nq S+1 . 

The following lemma collects auxiliary results which 
will become helpful later. 

Lemma 2 (Auxiliary results). 

1. LetveV. Then 

(v,s- V x v) L 2 (Q) > 3 f s f r _ (S) v 2 s ■ n(x)dxds. If 
furthermore v € Vo, then (v,s • ^ x v)l 2 (Q) — 0- 

2. Forv eH l -°(Q.), it holds \\s- V x v\\ < Vd\\v\\ H i, 0(Q) . 
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Proof. 1. A proof is given by (Manteuffel et al. 2000, 
Thm. 2.1). It uses the divergence theorem and 
exploits the fact that v|an_ = 0 for s • n(x) < 0 if 
v € Vo, where nix) is the outward unit normal on the 
boundary dD: 

(v,s- V x v) L t (Q) = \ I f V* • (sv 2 )dx ds 
2 Js Jd 

= - f f v 2 s ■ n(x)dx ds 
2 Js JdD 

= - f f v 2 s ■ n(x)dx ds 

2 Js Jr.(s) 

H — f f v 2 s ■ n(x)dxds. 
2 Js Jr+(s) 

As s • n > 0 in the second integral, we obtain the 
first assertion. If additionally v € Vo, then the first 
integral vanishes, and the second assertion follows. 
2. We again quote Manteuffel et al. (2000, Lemma 
4.1 (i)): 



lis • V*v| 



SiD Xi v I ds dx 



< d I I (siD Xi v) ds dx 

JdJs i=1 

d 

< d^2\\D xl vf < d\\vf Hh0(ay 



□ 



In order to establish well-posedness of the variational 
formulation (21), we prove continuity and coercivity of the 
bilinear form (19) and continuity of the linear form (20) in 
the following. 

Lemma 3 (Continuity of bilinear form). Let a, k, s e 
L°°(D) with \\a\\ L ^ {D) =: <x max , ||k|| £0 o (D) =: *r ma » 
Ikllioop) =: £ max> then there is a constant 0 < c c < oo 
such that for all u, v e V\ 

\a{u,v)\ < c c || m || i ||v||i. 

Proof. We proceed analogously to (Manteuffel et al. 
2000, Thm. 3.3). To begin with, we estimate for all u, v € V 

||Rv|| = ||e*rv + ss ■ W x v\\ < e m axf m axl|v|| + e ma x \\s ■ V x v\\ 
< max{e max /c max , 1} HvHx , (24) 
l|Aw|| < Kmaxllall + lis - V*w|| + 

"'max IIQl«ll 

}|||w|||. (25) 



Using the Cauchy-Schwarz inequality as well as esti- 
mates (24) and (25) it holds 



\a{u,v)\ < IHRvll ||A«||+2||v||_||«||_| 

<2(||Rv|| 2 +||v||2) 1/2 (||A M || 2 + || M || 2 _) 1/2 



< 2 maxjl, emax^max) max{ l.°Vr 



IlKlll IIVlll- 



□ 



Lemma 4 (Continuity of linear form). Given the assump- 
tions of Lemma 3 on k, a, s, and additionally f e L 2 {Q,), 
g : 9£2_ R with \\g\\_ < oo, there is a constant 
0 < C[ < oo such that for v € Vi it holds 

|/(v)| <c,||v||i. 
Proof. The proof is analogous to that of Lemma 3: 



1/2 



|/(v)| < | ||Rv|| |/||+2||v||_ \g\_\ 

<2(||Rv|| 2 +||v|| 2 _) 1 / 2 (| / | 2 +W 2 _) 
< 2max{l,e m axKmax}(||/|| + |^||_)l|v|li. 



□ 



Next, we show coercivity of the bilinear form. For ease 
of exposition we shall assume £ and k to be constant on 
the physical domain. Coercivity can also be obtained for 
non-constant coefficients (see Widmer 2009, Thm. 2.2, for 
an example). Coercivity of the SUPG variational formula- 
tion for the RTP has also been proved by (Avila et al. 2011, 
Lemma 2), although in a different norm. Previously, we 
had proved coercivity of the T-stabilized variational for- 
mulation without the boundary form b(; ■) (Grella 2013, 
Lemma 4.1), here we include this boundary form in the 
formulation, which will motivate the choice of the stabi- 
lization parameter e. 

Lemma 5 (Coercivity of bilinear form). Let k, b be posi- 
tive functions which are constant on the physical domain 
D. Assume mm x(z 0 o =: <T min > 0 an d er max = || er || L ao (D) , 



and additionally that 



2 2 
^max < * Ka min> S < 



(26) 



Then the bilinear form «(•, •) from (19) is coercive on Vi x 
Vi: there is a constant c e > 0 such that for all v e V\ it 
holds 

a(v,v) > c e \\v\\ 2 . 
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Proof. For an overview of the involved terms we split the 
bilinear form into separate inner products: 

a(v, v) = (skv + ss ■ W x v, Av) + 2b(v, v) 

= (skv + ss ■ W x v, s ■ W x v + kv + Qv) + 2b(v, v) 
= (skv, s ■ V x v) + (skv, kv) + (skv, Qv) 
+ (ss ■ V x v, s ■ V x v) + (ss ■ V x v, kv) 
+ (ss ■ W x v, Qv) + 2b(v, v) (27) 

As we assumed s and k to be constant, we can factor 
these coefficients out of the inner products. We begin by 
analyzing the sum of first and fifth inner product. 
Applying statement 1 of Lemma 2 yields 

2sk (v,s-V x v) > -sk\\v\\ 2 _. 
Together with the boundary term, we obtain 

(skv,s - W x v) + (ss - W x v,Kv)+2b(v,v) > (2-sk) \\v\\ 2 _. 
The second inner product is bounded from below by 

(skv,kv) > sk 2 \\v\\ 2 . 

To estimate the third inner product, we use prop- 
erty (11) of the scattering operator: 

sk (v, Qv) > sk \\Qv\\ 2 > SKal in \\Q lV f . 

The fourth inner product in Eq. (27) is 

(ss - V x v,s- V x v) = s \\s- V x v|| 2 . 

For the sixth inner product we apply Cauchy-Schwarz 
inequality and Young's inequality with a parameter 9 > 0: 

(£5 • V x v, Qv) > -ea max ||s • V*v|| ||Qiv|| 

'9 1 

> -Wmax| " l|S-V^|| 2 + — 

Combining all estimates yields the result: 



a(v,v) >£K 2 ||v|| 2 + £ 



(! 

ield 



IIQivll 



lis • V*v|| 



mm 20 

2 



IQivll 



+ (2-sk) ||v 

2 



> mm { sk 



> £ I J- 2 ^rnax I > 



S K<J, 



lnin " ^cr nax ),2-SK\\\v\\ l . 



By eliminating 9 we obtain the condition c r ^ lax < 
4/ctT min- ^ ne condition s < 2/k results from the last of the 
terms over which the minimum is taken. □ 



The previous condition on the stabilization parameter 
leads to a choice of e = 1/k. Well-posedness of the 
variational formulation now follows directly. 



Theorem 6 (Existence and uniqueness of solution to 
variational formulation). Provided that f e L 2 (Q.) and 
\g\_ < oo there exists a unique solution u e V\ to the 
variational formulation (21). 

Proof Since ( Vi , II - II i) is a Hilbert space and Lem- 
mata 3-5 guarantee continuity of the augmented SUPG 
bilinear form and linear form as well as coercivity of the 
bilinear form, the Lax-Milgram theorem (Brenner and 
Scott 2008, Thm. 2.7.7) ensures existence and uniqueness 
of the solution to (21). □ 

Discretization 

For the discretization of the variational problem (21), 
we restrict the space Vi in the variational formulation 
(21) to tensor products of hierarchic, finite dimensional 
approximation spaces over the component domains D 
and S. 

Full tensor discretization 

In the standard full tensor approximation, we choose a full 
tensor product space V L,N to approximate Vi: 



Vl w yL,N ._ yL 0 yN_ 



(28) 



As H lfi (Q) = H l (D) <g> L 2 (S) is a dense subspace of Vi, 
we define the family of physical approximation spaces as 

vjP := S 0 ' 1 (D, To D )cH l (D), l D = l,...,L, (29) 

the spaces of continuous, piecewise linear functions on a 
dyadically refined mesh Tjf over D. Here, the parameter 
Id stands for the physical resolution. It is related to the 
mesh width h in 7^° by h = O (2 _,D ). With respect to 

the resolution Id = 0, . . . ,L, the spaces V l D D form a nested 
sequence 

ygcv£c...cvgctf 1 (D). 

Let Md := dim denote the number of degrees of 
freedom for the FE space Vq in the physical domain D. 
Then 



M D = 0 (2 dL ^j 



(30) 



with the dimension d of the physical domain. The exact 
number will depend on the geometry of the domain. For a 
square or cube D =[ 0, l] d , respectively, we obtain 



M D = {2 L + l) d . 



(31) 



In the angular domain, we distinguish between the Pm- 
method and the SjV-method. 
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SAr-method Here, the family of approximation spaces is 
given by 

V l s s := S" 1 ' 0 (s, if) C L 2 (S), l s = l,...,N, 

(32) 

the spaces of piecewise constant functions on a dyadically 
refined mesh Tg S . As the physical spaces, these spaces 
exhibit a nested structure. The angular resolution N and 
the dimension of V$ are related by 



M s := dim Vf = O (V^) . 



(33) 



^Ar-method To define the angular approximation 
spaces of the PA^-method, we first introduce the spaces of 
spectral functions of the tig-sphere, 



N 



spanjY^ :n = 0,...,N; 

m = 1, . . . , m nids | CL 2 (S), 



(34) 



where Y$jp are the spherical harmonics of the d^-sphere, 
and m ni d s is the largest value of the secondary index m 
depending on the value of the primary index n and the 
dimension. These spaces offer an inherent nested struc- 
ture. To obtain the same relation (33) between resolution 
level and degrees of freedom as in the Sw-method, we 
connect the resolution level N and N by N = 2 N — 1. 
Then, the angular approximation spaces are 



V 1 /:- 



T> d S 

" 2*5-1' 



l s = h...,N, 



(35) 



and relation (33) also holds here. Up to the index rela- 
beling and the additional boundary form, we obtain the 
spherical harmonics method already analyzed by (Grella 
and Schwab 2011a). 

In both methods, the full tensor approximation space 
consequently has the dimension 

M L , N := dim V LM = M D ■ M s = O (2 dL+d s N \ . (36) 

The full tensor approximate solution can be expressed 
by means of a physical basis of Vq and an 

angular basis {/5j/}fcff of Vg as 

M D M s 

u l ,n(x, s) := Uijai(x)Pj(s) (37) 

with solution coefficients e R. The discrete variational 
formulation finally reads: Find uijq € V L,N such that 

a(u L , N ,v LM ) = l(v LM ) Vv L , N € V L,N , (38) 

with the bilinear form «(-, •) from (19) and the linear form 
/(•) from (20). As V L ' N is a subspace of Vi well-posedness 
ensured by Thm. 6 for the continuous problem follows 
also for this discrete problem. 



By choosing a subset of H (D) <g> L 2 (S) as trial space 
we effectively assume a slightly higher regularity on the 
solution than what is guaranteed by the definition (12) of 
Vi. For instance, solutions with line discontinuities due 
to the transport of discontinuous boundary data into the 
domain are not included in V L ' N . However, since V L,N is 
dense in Vi, even discontinuous solutions will be approx- 
imated with increasing resolution. Furthermore, in order 
to leverage the advantages of a sparse tensor approxima- 
tion, a higher regularity of the solution will be required in 
any case. 

Equivalence of collocation DOM and phase space Galerkin 
DOM with quadrature 

Ordinarily the discrete ordinates method is presented as 
a collocation method in angle: Fixed directions s ; € S, 
j = 1, . . . , Ms, are inserted into the RTE (la), and for each 
direction, the intensity Uj{x) := u{x,Sj) € Vk is sought 
as the solution to a purely spatial PDE. In these PDEs, the 
scattering integral is replaced by a quadrature rule 



, M S 

/ '3>(Sj,s')u(x,s')ds' % w m ${Sj,s m )u, 

^ m=\ 



(39) 



with weights w m > 0. By applying a Galerkin ansatz with 
stabilization in the physical domain to the PDEs, a system 
of coupled variational formulations 



I R ; v, TjUj + auj - ^2 w m s m )u m J 

\ m=l / 



L 2 (D) 



+ 2(v, \sj-n\uj) LHr _ isj)) 

= M)lHd) + 2 ( v > I*/ ' "tiW-W Vv 6 V D 

(40) 

results with directional stabilization and transport 
operators 

Ry:=R| s=s/ , T ; :=T| S=S ., j=l,...,M s . (41) 

In the phase space Galerkin approach, variational for- 
mulation (38) is discretized further by substituting the 
angular quadrature rule (39) for all angular integrals so 
that the bilinear form (19) is approximated by 

M S 1 

a(u, v) « a(u, v) = ^ Wj I R ; v/, TjUj + auj 
j=i \ 

M S 

~ ^2 W m ${Sj,S m )u m 

m=l / /J!(D) 

M S 

+ 2 J £ i W j (vj,\s r n\uj) LHris/)) . 

7=1 
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Let the linear functional /(•) from (20) be approximated 
by a functional /(•) with angular quadrature correspond- 
ingly, then the directional solutions Uj are determined 
from the variational formulation with angular quadrature 

a(u, v) = 7(v) Vv € V L ' N , j = 1, . . . , M s . (42) 

Since this formulation has to hold for all v e V L,N , 
it follows that for test functions which vanish at every 
angular quadrature node s,-, i = 1, . . . , Ms except one Sj, 
formulation (42) can be reduced to the variational for- 
mulation (40) from the collocation discretization. This 
condition on the test functions is satisfied e. g. for a basis 
of the test space of characteristic functions on the angu- 
lar mesh if each mesh cell contains exactly one angular 
quadrature node. With such a one-point quadrature rule 
and characteristic basis functions of Vg , the phase space 
Galerkin DOM is therefore equivalent to the collocation 
DOM after discretization. 

Sparse tensor discretization 

The full tensor approach presented before shows the typ- 
ical complexity for full tensor approximations: The num- 
ber of degrees of freedoms increases exponentially with 
the dimension and the resolution levels in a dyadically 
refined scheme. 

A way to counter this exponential increase is found 
in sparse tensorization. Using the same approximation 
spaces on the component domains Vp and V 1 / as for 
the full tensor approximation we define a sparse tensor 
approximation space V L ' by 

Vl « v L - N := V D ® V S> (43) 

o</(Ws)<£ 

where the sparsity profile f : [0,...,L] x {0, . . . ,N] — »■ K 
determines which tensor product subspaces Vp <g) V 1 / are 
to be included in the approximation. The sparsity profile 
usually depends on N as well. Here, we employ a linear 
profile 

f(l D ,ls) = l D +Ll S /N, (44) 

which is normally chosen if the component complexities 
Mo and Ms depend on the resolution parameters L and 
N in the same way and identical order of approxima- 
tion is sought over both component domains (cf. Zenger 
1991; Bungartz and Griebel 2004; Griebel and Harbrecht 
2013a). 

If direct sum decompositions of the component approx- 
imation spaces V D D and vlf into detail spaces Wp and 
W l /,i.e. 

V l D D = V l D D - 1 ®W D D , l D = l,...,L 



are available (correspondingly in the angular domain), 
then the sparse tensor approximation space V L,N can also 
be written as 

yL,N = J2 W£®wjf. (45) 
0<f(l D ,l s )<L 

By choosing hierarchical bases for V D D and Vg 8 , each 
degree of freedom tig can directly be associated with a ten- 
sor product detail space W D D <g> W s s . The sparse solution 
is then given by 

UL.N = E U[ Di [ s , 

0</(to,/ s )<£ 

dim VPj° dim W l £ 

«w s = E E ¥(Wf(^<^ s 

(=1 ;=1 

Thus, the sparse discrete variational problem reads: Find 
ul,n e V LM such that 

a{u L , N ,v LiN ) = 1(v l ,n) Vvl,n G V lm . (46) 

The dimension of the sparse tensor product space V L,N 
depends on the sparsity profile /(fo, Is)- For a linear spar- 
sity profile as in (44), the following complexity estimate is 
known (e. g. Bungartz and Griebel 2004, Lemma 3.6), or 
Griebel and Harbrecht (2013a, Thm. 4.1)). 

Lemma 7. Assuming the dimensions of the detail spaces 
W l g and W 1 / scale as dim(wf') < c{l dik with constants 
d > 0 and dimensions dj, i = D, S, with do = d, and given 
a linear sparsity profile f(lr>, Is) as in (44), the dimension 
of the sparse tensor product approximation space V L,N as 
defined by (45) is 

M Li n < L 9 2 max[dL ' dsN] < (hgM D ) d max{M D ,M 5 }, (47) 

where 0 = 1 if dL = d$N and 9 = 0 otherwise. Relation 
"<" defines an order up to constants with respect to the 
relevant scaling parameters L, N: a < b iff a < Cb with 
constant C independent ofL and N. 

Error analysis 

In this section, we shall show that the convergence rates 
of the full tensor and sparse tensor Galerkin methods dif- 
fer only by a logarithmic factor in the degrees of freedom, 
provided that somewhat stronger regularity requirements 
are met for the exact solution. 

The analysis will proceed along the usual fashion, cp. 
(Bungartz and Griebel 2004). We define the Galerkin pro- 
jector P L - N : Vi V L ' N into the full tensor product 
approximation space 

a (V L ' N u, v) = a(u, v) Vv € V L,N . (48) 
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Letting L —*■ oo (N oo) the fact that the sub- 
spaces are closed and dense implies that in the respec- 
tive limits we obtain semidiscrete Galerkin projectors 
P£ := limz^oo P L ' N (pjj := lim^^ V L > N ) on the physi- 
cal (angular) domain, as the Galerkin projector is stable in 
the || -||i -norm: 

Lemma 8 (Stability of the Galerkin projector). Let v e Vi. 
Then there is a constant cp > 0 independent ofL and N so 
that 

\v LM v\i <<*IMIi. 

Proof. With continuity (Lemma 3) of the bilinear form 
we obtain 

\a(P L - N v,v LrN ) | = \a(v,v L>N )\ 

< C c ||v||l || V iiA f || ! Vvi,tv G V L,N . 

Since this holds for all vl,n e V L,N , we can set 
v LtN = P L ' N v and exploit coercivity of the bilinear 
form (Lemma 5): 

Ce ||P L '^||J< ^(P^V.P^V)! 

= l«(v,P i ' N v)|< Cc ||v|| 1 IP^HIi. 

If p^v / Owe obtain the result with cp = c c /c e . □ 

Error estimates on the physical domain 

To begin with, we require some approximation results in 
the// 1 (£>) 

-norm on the physical domain. With a Clement- 
type quasi-interpolation operator Pj (Scott and Zhang, 
1990, Thm. 4.1 and Cor. 4.1) we obtain 

Lemma 9 (Approximation of quasi-interpolation). For 
polyhedral D C K rf and a shape-regular triangulation 
on D with mesh width h = 2~ L , the quasi-interpolation 
P[v of a function v e H S+1 (D), s e [0,1], to the space 
Vp = S 0 ' 1 (D, 7q) of piecewise affine functions on T£ 
satisfies the error estimate 

||v-Pfv|| H i (£)) < c H 2- sL \\v\\ H s + i (D) , 

where ch > 0 is a constant independent ofL. 

Lemma 10 (Stability of quasi-interpolation). Under the 
assumptions of Lemma 9, quasi-interpolation is H l -stable, 
i. e. there exists a constant eg > 0 independent ofL such 
that for all v e H 1 (D) it holds 

l|Pfv|| H i (D) < cb||v|| h i (£)) . 

Next we derive an error estimate for the Galerkin 
approximation on the physical domain. At this point, the 
approximation is semidiscrete. 

Lemma 11 (Error estimate for Galerkin projection on 
physical domain). Let u e Z/ s+1 '°(fi), s e {0,1}, be the 



exact solution to problem (21) and Ui := Pj^u e ® 
L 2 (S) the Galerkin projected solution to 

a(u L , Vi) = l(v L ) Vvi e Vj§ <8> L 2 (S) (49) 

with «(•, •) from (19) and /(•) from (20). Then, there is a 
constant c p > 0 independent ofL such that 

II" -Willi < c p 2~ sL \\u\\ H s+ifi (n) . 

Proof. The proof is standard, and is based on coerciv- 
ity and Galerkin orthogonality. We proceed analogous to 
(Avila et al. 2011, Lemma 3 and Theorem 1). After insert- 
ing the quasi-interpolated solution ui := (P[(giId l s)M with 
P[ from Lemma 9 the triangle inequality permits us to 
write 

llw - Willi < || M - "if! + II" 1 ~ Mi lli • ( 5 °) 

For the first part, we use the fact that there is a constant 
c„ > 0 for all v e H l (D) <g> L 2 (S) such that 

IMIi < c n \\v\\ H i,o (a) . 

Thus, we can apply Lemma 9: 

||w — < c n \ u — ui | w i,o (n ) < c„c H 2~ sL \\u\\ H s+i,0( n) . 

For the second part in (50), we use coercivity of the bilin- 
ear form, then in a second step Galerkin orthogonality, 
and finally continuity of the bilinear form to write 

\\u L - Mi|| j < c~ l a{u L - u L ,u L - u L ) 

< c~ l a(u - ul, ul - ul) 

< c c c~ l J u - ul 1 1 I u L -ul\ x 

< C c C~ l C„ \u - U L \\ Hh0(a) \\u L - Wi|| j , 

and therefore with Lemma 9 

||wi - Wifi < c c c~ l c n c H 2~ sL \\u\\ H s+i,o (Q) . 
By inserting into (50) we arrive at the result 

II W - Willi < C n CH (1 + C c C~ l ) 2~ sL \\u\\ H s+l,o (Q) . 

□ 

Error estimates on the angular domain 

On the angular domain, the considerations in the follow- 
ing require an approximation result for Z, 2 -projections. 

Lemma 12. For functions v e H t (S), t e {0,1}, the L 2 - 
projection to the space satisfies the error estimate 

\\ v -^H S A L H S) - Cl2 ~ tNMH '^' (51) 
where the constant c/ > 0 is independent ofN. 

This result can be obtained for approximation by spec- 
tral functions as in the spherical harmonics method (in 
which case t > 0 is arbitrary), for instance, as well as for 
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approximation by piecewise constants as in the discrete 
ordinates method (in which case 0 < t < 1). It allows 
the derivation of the same approximation rate for the 
semidiscrete Galerkin projection on the angular domain. 

Lemma 13 (Error estimate for angular Galerkin projec- 
tion). Let u 6 H l,t (£l), t € {0, 1}, be the exact solution to 
problem (21) and um '■= P^w e H l (D) <g> the Galerkin 
projected solution with angular part from the subspace Vg 
ofL 2 (S). Then there is a constant c a > 0 independent ofN 
such that 

\\u-u N \\i <c a N- 1 \\u\\ H i, t(Q) . 

Proof. The proof proceeds analogously to the one of 
Lemma 11 while substituting the L 2 -projected solution 
with Lemma 12 for the quasi-interpolated solution, the 
details are therefore omitted here. □ 

Error estimate for the full tensor phase space Galerkin method 

The following theorem gives an error estimate for the full 
tensor approximation. 

Theorem 14 (Error estimate full tensor Galerkin 
method). The full tensor Galerkin approximation 
u LM = P L ' N u of a solution u e H s+1 >°(£i) n H l '\Q.), 
s € {0,1}, t € {0,1}, to the variational problem (21) 
satisfies the asymptotic error estimate 

\\u - u LtN \ l < 2~ sL ||k|| ws+ i,o (£2) + 2~ tN \\u\\ H u (n) , 

(52) 

with relation "<"as in Lemma 7. 

Proof. By Cea's Lemma (Brenner and Scott 2008, Thm. 
2.8.1) the Galerkin approximation is quasi-optimal in 
V L,N , its error can therefore be bounded (up to con- 
stants) by the error of any other approximation to u in 
V L,N , for example the quasi-interpolated andi 2 -projected 
approximation P[ (g> P^ 2 w: 

||M-P i * Af w| 1 < ||w-Pf (giP^w^ < ||w-Pf (gildwlj 
+ ||Pf <g)Idw-Pf (giP^wlj 
<2- sL \\u\\ H s +1 ,o (a) 
+ I (Id - Id <g> Pf 2 )Pi <E) Idw || x 

< 2~ sL \\u\\ H s +1 ,o (a) 

+ 2~ tN |P[ ®ldu\\ HU{n) 

< 2~ sL ||m|| hs +i,o (S2) + 2~ tN \\u\\ H U(Q) ■ 

Here, we used the approximation properties of the 
quasi-interpolant from Lemma 9 and of the angular L 2 - 



projection from Lemma 12. The last step is a conse- 
quence of the H 1 -stability asserted in Lemma 10 of the 
quasi-interpolation. □ 

Error estimate for the sparse tensor phase space Galerkin 
method 

After the full tensor approximation properties, we con- 
sider the convergence properties of a direct sparse tensor 
approximation on the sparse tensor product space V L,N as 
defined in (45). 

In analogy to the full tensor Galerkin projector P ^ , we 
can define a sparse tensor Galerkin projector P L,N by the 
orthogonality relation 

a (p l - m u, v) = a(u, v) Vv € V L ' N . 

The error of the sparse tensor solution 
is estimated in the following theorem (see also Widmer 
2009, Thm. 2.6) and (Griebel and Harbrecht 2013a, 
Thms. 4.3 and 7.1)). 

Theorem 15 (Error estimate of direct sparse tensor solu- 
tion). Let the linear sparsity profile as in (44) be given. 
Assume further that L and N vary such that —sL + tN = 
f = const, then the direct sparse tensor approximation 
ul,n °f a function u € // s+1,t (£2), s,t <= {0, 1}, satisfies the 
error estimate 

\u ~ Mi^lj < L (2- sL +2~ tN ) |M|„ s+ « ( r2) , 

where relation "<" is defined as in Lemma 7. 

Proof. We follow the proof of Thm. 2.6 by (Widmer 
2009). First we introduce so-called difference projectors 
a Id ._ p fo _ pfo-i and A ls ._ p fe _ pjs-i as the dif . 

ference between projections to two consecutive resolution 
levels with the convention Pj" 1 = 0 = P~ 2 \ They project 

onto the detail spaces W l § and W l £ , respectively. 

With these difference projectors, a sparse quasi- 
interpolated andi 2 -projected approximation iti^r G V L,N 
to u can be expressed as 

L lT X Co) 

UL,N =Y1 J2 A I° 0 ^l? 11 ' 
l D =0 l s =0 

where 1™ x (Id) is the largest feasible angular resolution 
index which results from solving /(t>, Is) < L with 
respect to lg. 

Now we exploit quasi-optimality of the Galerkin 
approximation on the sparse tensor product space to 
replace the Galerkin approximation error by the error of 
the quasi-interpolated and i 2 -projected approximation. 
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Additionally applying the norm estimate ||v||i < 
yields 



| U - Ul,N II j 



< 



« - E E A i D ® a1 3 u 

fo=0 ls=0 



« 1 .°(Q) 

(53) 



The error is split into two terms: 



|| « "1,^11^1,0(0) 



E E a;-®a^ 

fo=o/ s =fs" x (fo)+i 



+ 



oo oo 



E 

; D =i+i i s =o 



(54) 

The second term on the right hand side can be estimated 
by Lemma 9: 

// = || (Id - Pf) ® Idu\\ Hh0(Q) < c H 2~ sL \\u\\ H s+i,o (Q) . 

(55) 

This term will not contribute to the asymptotic terms. 
The first term on the right hand side of (54) is split up 
further: 



^(pfe_ P ; D -i) 0 ( Id _ P ™), 

'd=0 

£ (P<» - Id + Id - Pj-" 1 ) ® (id - P ! p M ) , 
'd=0 

£(l("-*W H -'r>U, 



H li0 (Q) 



b=0 
+ 



(56) 



Both norms on the right hand side of (56) can be 
estimated by Lemma 9 and Lemma 12: 

< c H T slD II Id ® (id - pJ X(/d) ) u\ 

<c H c{T dD - tl ™ {lD) \\u\\ H s + u (Q) . 
Inserting back into (56) yields 



I<2c HCl \\u\\ H s + u (a) J22- slD - tl ^ (lD) . 
i D =o 



(57) 



The task is now to estimate the series. Using the 
assumption £ = — s + tN/L: 

L L 

2- sl o-tN/L(L-l D ) _ 2~tN 2(- s + tN / L '> l D 

l D =0 l D =0 

L 



-tN 



E 2 ffo . 



(58) 



We estimate the sum on the right hand side of (58) by its 
largest summand. Two cases can be distinguished here: 

1. If f < 0, the largest summand occurs for Id = 0: 

L 

2 -tN J2 2 (lD < L2~ tN . 

l D =0 

2. If £ > 0, the largest summand occurs for Id = L: 

L 

2~ tN 2^' D < 2~ tN L2~ sL+tN = L2~ sL 
l D =o 

In summary, we may write 

L 



2 - s 'o- tf s ax (fc) < 12' 

i D =o 



■sL—tN 



By combining this estimate with relations (53) to (57), 
we finally arrive at 

< L2~ sL - tN I 



\u- u LtN \\ Hh0(n) 



\U\\ H s+l,t (n) . 



□ 



In conclusion, we find that the convergence rate of 
0(2~ sL ~ tN ) of the full tensor approximation is maintained 
up to an additional factor L, which by Md = O (2 dL ) is log- 
arithmic in the number of degrees of freedom. This result 
in conjunction with the greatly reduced complexity of the 
sparse tensor method (Lemma 7) shows its superior effi- 
ciency provided that the function u to be approximated is 
at least in H S+U (Q.), with s, t e {0, 1}. 

Numerical experiments 

Algorithms 

For the numerical experiments we compute a sparse ten- 
sor solution with the help of the combination technique. 
The sparse solution is constructed according to the for- 
mula 



UL.N ■ 



E ( Ue »J 



from a number of solutions Ui Dl i s € V 1d,1s to the full 
tensor discrete variational formulation (38) of reduced 
physical resolution Id and angular resolution 1$. 
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Clearly u m is in the space V L - N = £[ D=0 v to ^s' <fD \ 
which is identical to the sparse tensor approximation 
space from (43). However, in general the combination 
approximation differs from a direct sparse approxima- 
tion uljj (see also Grella 2013, Sec. 2.3.1). Due to 
the quasi-optimality of the direct sparse solution as 
an approximation in V L,N , the error of the combina- 
tion approximation can serve as an upper bound (up to 
factors) for the error || u — ui,u \ x of the direct sparse 
approximation. 

Note that the convergence of the combination technique 
for the radiative transfer problem has not been shown for- 
mally yet. A recent proof for elliptic operators by (Griebel 
and Harbrecht 2013b) would be applicable under certain 
stability assumptions on the semidiscrete Galerkin pro- 
jectors (for details we refer to (Grella 2013, Sec. 5.3.7)). 
However, the use of the combination technique approx- 
imation has practical advantages over the direct sparse 
approximation. First, to construct the subproblem solu- 
tions of lower resolution, an existing full tensor solver 
with standard nonhierarchical FEM bases can be reused, 
no direct sparse solver needs to be implemented. Second, 
the splitting into subproblems entails a natural level for 
parallelism in the algorithm, which can still be combined 
with parallel solution procedures at the level of each sub- 
problem (an implementation is described in (Grella 2013, 
Chap. 7)). 

Each of the full tensor subproblems is solved by a 
phase space Galerkin finite element method with non- 
hierarchical affine hat functions as physical basis and 
piecewise constants as angular basis. In the experiment of 
Sec. 'Experiment 2', the midpoint rule is used for angular 
quadrature which corresponds to the %-method. How- 
ever, in situations where ray effects (Lathrop 1968) pollute 
the results, adaptive quadrature may help (Stone 2007). 
As a simple adaptive rule we link the number of quadra- 
ture points n q per dimension and per mesh element to 
the resolution levels Id, Is of the subproblem by n q = 
max{/fl//s, 1} in the experiment of Sec. 'Experiment 1! 
Even though the overall computational effort is then not 
bounded by Lemma 7, the total number of degrees of free- 
dom still is. As the iterative, approximate solution of the 
linear system constitutes the most time consuming part, 
the sparse tensor method is, in practice, more efficient 
than the full tensor method. 

Quantities of interest 

In applications, the radiative intensity is often coupled to 
other modes of energy transport via the net emission (e. g. 
Larsen et al. 2002, Eq. (1.1a)). The net emission can be 
computed in turn from the incident radiation 

G(x) = j u(x,s)ds. (59) 



For this reason, we choose the incident radiation as a 
lower-dimensional variable to visualize results and to ana- 
lyze errors. The relative L 2 - or // 1 -error of the incident 
radiation is given by 

err(G L , N ) x = \\G - G LtN \\ x /\\G\\ x , X = L 2 (D),H l (D). 
Numerical experiments 

All experiments are set on the domains D = [0, S = 
S ds , with d = ds + 1. We solve the RTP with isotropic 
scattering <I>(s,s') = 1/|5| and zero inflow boundary 
conditions g = 0. 

Experiment 1 

We search the solution to the Gaussian blackbody radia- 
tion 

I b (x) = 2exp (-3 2 (* - c) 2 ) , c = (0.5,0.5) T , 

with absorption and scattering coefficient k = a = 1. 

The H 1 -error of the incident radiation indeed con- 
verges faster in the sparse approximation than the full 
approximation (Figure 1). Note that the L 2 -error of the 
sparse approximation can be larger than the error of the 
full approximation because the sparsity profile /(Id, Is) 
has been optimized for essentially undeteriorated conver- 
gence in the ||-||i-norm of the error in the radiative inten- 
sity, which is more closely represented by the H 1 -error 
than the L 2 -error of the incident radiation. 

Experiment 2 

A blackbody radiation Ib(x,s) corresponding to the exact 
solution 

3 3 
u(x,s) = — (1 + (s ■ s' f ) Y[(-4Xi(Xi - 1)), 

i=i 

with fixed s' = (l/\/3> l/\/3> l/\/3) T is inserted into the 
right hand side functional in (38) (Grella 2013, Sec. 8.2, 
Exp. 1). The absorption is set to k = 1, the scattering 
coefficient to a = 0.5. 

For this experiment we employed a discrete ordinates 
solver in which the angular resolution N' is related to the 
angular degrees of freedom by Mg = (N' + l) 2 so that 
N «s Llog 2 (AT' + 1)J, where N is the angular resolution 
used otherwise in this paper. 

Figure 2 shows the superior efficiency of the sparse 
approach with respect to number of degrees of freedom 
vs. achieved error. The convergence rates indicate that the 
curse of dimensionality is mitigated by the sparse discrete 
ordinates method. 

For a comparison to other sparse tensor approaches 
we refer to the numerical experiment section of (Grella 
2013), which features a sparse tensor spherical harmon- 
ics approximation and a sparse collocation discrete ordi- 
nates method realized via the combination technique. We 
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2D Gaussian, L = {1 ,2,3,4,5}, N = {2,3,4,5,5} 

10° e — — — 




Total #DoFs 



Figure 1 Experiment 1: Convergence in incident radiation with full and sparse phase space Galerkin approximation. Reference resolution 
was i re f = 6/W re f = 6. Reference slopes provided as visual aids only. Even with the lowest order sparse tensor phase space Galerkin discretization, 
the savings in DoFs to reach engineering accuracy of 1 96— 1 096 in the H ] error is about an order of magnitude. 



observed that the approach presented here performs sim- 
ilarly to the sparse collocation DOM combination tech- 
nique as the methods are similar from the point of view 
of implementation, even though their theoretical deriva- 
tion is different. The presented approach is somewhat less 
susceptible to ray effects at the expense of slightly longer 



computational times as the angular quadrature is adapted 
to the resolution of the angular mesh. The spherical har- 
monics method is most effective for solutions with highly 
regular angular part because of its regularity requirements 
for spectral convergence. In general, at the same resolu- 
tion levels L and N, the combination technique approach 




Total #DoFs 



Figure 2 Experiment 2: Convergence in incident radiation with full and sparse DOM. Reference resolution was L le f = 4. Angular resolution N' 
corresponds to N ~ {1,2, 3,4). Reference slopes provided as visual aids only. The savings in DoFs to reach engineering accuracy of 1 96—1 096 are 
about two orders of magnitude. 
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realizes approximately the same error as the direct sparse 
approach, while the number of degrees of freedom in the 
combination technique is larger than in the direct sparse 
approach because the approximation spaces of different 
subproblems in the combination technique overlap in the 
degrees of freedom. It is therefore slightly less efficient 
than the direct sparse approach, but considerably more 
efficient than the full tensor approach and advantageous 
in practice due to faster and simpler implementation and 
parallelization. 

Conclusion 

We have shown a direct sparse tensor phase space 
Galerkin approximation of the radiative intensity in the 
stationary monochromatic radiative transfer problem can 
be computed with only OilogMoiMo + Ms)) degrees 
of freedom as opposed to O(MoMs) degrees of freedom 
for a standard full tensor approximation. Here, is the 
number of physical degrees of freedom and M$ the num- 
ber of angular degrees of freedom. At the same time, the 
error of the sparse approximation in the ||-||i-norm still 
decreases essentially as the error of the full approximation, 

namely with the order O (logM D + M~ t/ds ^ 

as compared to O (M- s/d +M- t/ds \ in the full tensor 

approximation. The parameters s,t e {0,1} indicate the 
regularity of the exact solution which is required to be in 
the space of mixed smoothness H s+1,t (D x S) to achieve 
the sparse convergence rate, whereas H S+1 '°(D x S) Pi 
H l,t (D x S) is sufficient in the full tensor approximation. 

To simplify implementation, we realized the sparse ten- 
sor approximation algorithmically via the combination 
technique. Together with suitable quadrature rules, we 
demonstrated in numerical experiments that this sparse 
tensor combination approximation retains the analyzed 
theoretical advantages of the direct sparse tensor method 
while allowing for straightforward parallelization also at 
the level of subproblems. 

The proposed specialization of the phase space Galerkin 
framework investigated here has the advantage that both 
discrete ordinates and spherical harmonics method can 
be derived from it so that the sparse tensorization 
benefits hold for the sparse variants of both methods 
alike. 

Therefore, for problems whose solutions exhibit so- 
called mixed regularity, the sparse tensor product 
phase space Galerkin approximations realize a significant 
increase in efficiency, i.e. achievable error per number 
of degrees of freedom. Even in applications where high 
numerical accuracy is the main objective, a sparse ten- 
sor product approximation might be of value as an initial 
value for an iterative solver or in a problem-adapted pre- 
conditioning scheme. 
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